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ABSTRACT 

We develop an analytic theory to describe spiral density waves propagating in a shear- 
ing disc in the weakly nonlinear regime. Such waves are generically found to be excited 
in simulations of turbulent accretion disks, in particular if said turbulence arises from 
the magneto-rotational instability (MRI) . We derive a modified Burgers equation gov- 
erning their dynamics, which includes the effects of nonlinear steepening, dispersion, 
and a bulk viscosity to support shocks. We solve this equation approximately to ob- 
tain nonlinear sawtooth solutions that are asymptotically valid at late times. In this 
limit, the presence of shocks is found to cause the wave amplitude to decrease with 
time as t~ 2 . The validity of the analytic description is confirmed by direct numerical 
solution of the full nonlinear equations of motion. The asymptotic forms of the wave 
profiles of the state variables are also found to occur in MRI simulations indicating 
that dissipation due to shocks plays a significant role apart from any effects arising 
from direct coupling to the turbulence. 
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1 INTRODUCTION 

Accretion discs are common in astrophysics, occurring in 
close binary systems, active galactic nuclei and around pro- 
tostars where they provide the environment out of which 



planets form (see e.g. Papaloizou & Lin 1995 
paloizou 1996, for reviews 



Lin & Pa- 



Observational inferred accretion rates require enhanced 
angular momentum transport to take place. This is believed 
to be mediated by some form of turbulence, which can 
be regarded as providing an effective anomalous viscosity. 
The most likely source of accretion disk turbulence is the 



magneto-rotational instability (MRI, see Balbus & Hawley 



1991 19981. 



In both local and global simulations of the MRI with 
and without net magnetic flux, prolific density wave excita- 



tion has been found (e.g. Gardiner & Stone ( 2005[ ) in the 
local case, and Armitage (19981 in the global case). Such 



waves may be important in the contexts of quasi periodic 
oscillations. They may also be associated with stochastic 
migration of protoplanets ( |Nelson fc Pa paloizou 2004 Nel- 
|son|2005[ ) and with providing some residual angular momen- 
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turn transport in magnetically inactive dead zones [G"ammie| 
(1996 1. 



In view of the significance of these waves it is impor- 
tant to gain an understanding of the processes leading to 
their excitation and ultimate dissipation. In Heine mann fc| 



Papaloizou (2009a I we developed a WKBJ theory for the 



excitation of these density waves in the linear regime and 
found that vortensity fluctuations are responsible for their 
excitation during a short period of time, as they change from 
being leading to trailing, denoted as swing. In a second pa- 
per ( Heinemann fc Papaloizou|2009b I, we studied the wave 
excitation process directly as it occurs in fully non-linear 
three dimensional numerical MRI simulations. The results 
were found to be in good agreement with the WKBJ de- 
scription of the excitation process developed in [Heinemann| 
|fc Papaloizou] ( |2009a| ). 

In this paper we extend the work of our previous two 
papers to consider the behavior of the waves as they enter 
the nonlinear regime. This is significant because the even- 
tual development of weak shocks is expected and this leads 
to the dissipation of the waves. Their maximum amplitude 
can be expected to be determined by a balance between ex- 
citation and dissipation processes. The latter might occur 
either directly through wave phenomena such as shocks, or 
through interaction with the turbulence. Possibilities here 
include energy loses resulting from random secondary den- 



© 2011 RAS 



2 T. Heinemann and J. C. B. Papaloizou 



sity wave generation through interaction of the primary wave 
with turbulent eddies (e.g. Lighthill 1953 Howe 1971a|b 



Ffowcs Williams fc Howe|1973 l. Another phenomenon that 
could potentially play a role is distortion of the primary 
wave front on account of a variable wave propagation speed 
caused by the presence of dynamic turbulent eddies (e.g. 
Hesselink fc Sturtevant|1988[). Understanding these phenom- 



ena is notoriously difficult and their theoretical description 
is, at present, somewhat imprecise and associated with is- 
sues of interpretation. In the body of this paper we focus on 
weakly nonlinear wave theory without turbulence and de- 
fer consideration of wave-turbulence interactions until the 
discussion section. 

Our studies, in common with the majority of those 
that incorporate turbulence resulting from the MRI, are re- 
stricted to local isothermal disc models. Accordingly they 
do not incorporate the possibility of the channeling of wave 
energy towards the upper disc layers, that may occur for 
density waves with short wavelengths in the direction of 
the shear, in models with significant thermal stratification 



(Ogilvie & Lubow|1999 


i. 


As shown in 


Heinemann & Papaloizou 



bulent vortensity fluctuations arc able to excite pairs of 
counter-propagating spiral density waves when they swing 
from leading to trailing. Subsequently, as the excited waves 
propagate, their wavelengths in the direction of propagation 
decrease and the conservation of wave momentum results in 
an increasing amplitude, causing them to enter the nonlinear 
regime in which the formation of shocks occurs. Our aim is 
to describe the time asymptotic form of these excited waves 
and the way in which the associated shock dissipation leads 
to amplitude decay, and to relate these phenomena to what 
is seen in MRI simulations. 

We consider the nonlinear evolution of a single pair of 
counter-propagating spiral density waves on a homogeneous 
background. As we will demonstrate, this becomes an effec- 
tively one-dimensional problem when described in a shearing 
coordinate frame, reducing to ordinary one-dimensional gas 
dynamics in the far trailing regime. 

The plan of the paper is as follows. In Section [2] we de- 
scribe the shearing box model and present the basic equa- 
tions. In Section [2.1 1 we introduce shearing coordinates, giv- 
ing the form of the basic equations governing solutions that 
are functions only of the shearing coordinate in the direc- 
tion of the shear and time. We go on to express these using 
a Lagrangian formalism in Section [2.2| deriving the appro- 
priate forms of specific vorticity or vortensity conservation 
in Section 



2. I 



and nonlinear wave momentum conservation 



in Section 

In Section [3] we use the Lagrangian formulation as the 
basis for the development of an analytic description that is 
applicable to the weakly nonlinear regime. In this formula- 
tion, third order effects such as entropy production in shocks 
are neglected, making it possible to adopt an adiabatic or 
isothermal equation of state (e.g. |Yano|1996| |. In Section [3~T 
we derive a modified Burgers equation governing the unidi- 
rectional propagation of weakly nonlinear waves (for a sim- 
ilar approach to the problem of nonlinear planetary wakes 
Goodman & Rafikov 2001). This is found to provide the 



correct nonlinear development of wave profiles in which weak 
shocks arc ultimately present. Nonlinear sawtooth solutions, 
which are asymptotically valid at late times, are derived in 



Section |3.2| Corrections to the wave profile arising from dis- 
persive effects are derived in Section [3. 3| 

In Section [4] we go on to present numerical solutions 
of the full nonlinear equations governing a pair of counter- 
propagating spiral density waves undergoing swing excita- 
tion. At late times these are compared with analytic solu- 
tions derived in Section |3.2| The form of the wave profiles 
and the decay of the wave amplitudes with time are found 
to be in excellent agreement. Finally in Section [5] we dis- 
cuss our results, comparing them with what is seen in, and 
considering their consequences for, MRI simulations. 



MODEL BASIC EQUATIONS AND 
PRELIMINARIES 



We adopt the standard local shearing box of |Goldreich fc| 
Lynden-Bell ( 1965 1 and assume that the effect of magnetic 



fields on the density waves we consider is negligible. The 
basic equations are the conservation of mass and momentum 
in the form 



dp 

dt 

dv 
dt 



+ V • (pv) = 



■ v ■ Vv + 217 x v 



VP 

P 



2qQ, xe x + f. 



(1) 



(2) 



Here p, P, v, and / are the density, pressure, velocity, and 
viscous force per unit mass, respectively. We adopt a Carte- 
sian coordinate system (a;, y, z) with origin at the center of 
the box. The a;-direction points away from a putative cen- 
tral object, the y-direction is in the direction of rotation 
and the z-direction is parallel to the rotation axis. The unit 
vectors in these coordinate directions are e x , e v , and e z re- 
spectively. The angular velocity of the system is CI = Qe z . 
For a Keplerian disk the constant q — 3/2. 

In order to complete the description, we require an 
equation of state. This may be found by noting that the 
pressure P can be taken to be a function of the density p 
and the entropy per unit mass S. From the first and second 
laws of thermodynamics, the entropy satisfies the equation 



95 
dt 



+ v-VS 



(3) 



where T is the temperature and e is the net rate of heating 
per unit mass. When the latter is zero - as we will assume 
throughout this paper - then the motion is adiabatic, with 
S being conserved on a fluid element, and we may write S = 
So, where So is the fixed entropy per unit mass. We point 
out that heating by shocks will, in general, cause departures 
from an adiabatic equation of state, even in the absence of 
cooling. However, such effects enter only at third order in the 
wave amplitude (e.g. Yano|1996 1 and hence do not affect the 
dynamics of weakly nonlinear waves of the type considered 
in this paper. 

The basic state on which the waves we eventually con- 
sider propagate is such that P, p and S are constant and the 
velocity v = —qQxe y . Then for adiabatic motion So is con- 
stant and the relation P = P(p, So) leads to an effectively 
barotropic equation of state. We neglect vertical gravity and 
thermal stratification and look for waves that have no depen- 
dence on z. We note, however, that in the strictly isothermal 
case (P oc p), the velocity amplitudes associated with the 



© 2011 RAS, MNRAS 000,[T]{T3] 



Weakly nonlinear spiral density waves 3 



density waves we consider turn out not to depend on z even 



if the disk is vertically structured (Fromang & Papaloizou 



2007). But note that this situation could be modified for 



density waves with short wavelength in the direction of the 
shear if the disc model were to have significant thermal strat- 
ification ( |Ogilvie fc Lubow|| 1999 k 



2.1 Shearing coordinates 



In the theory of linear wave excitation developed in |Heine-| 
|mann fc Papaloizou ( 2009a b ) the coordinate dependence of 
the asymptotic form of the excited waves at late times is 
through an exponential factor 



exp[ifc a (y + qQtx) + ik x x^ , 



where k y is the constant azimuthal wavenumber and k x is 
the radial wavenumber at t = 0. This latter quantity can be 
removed by redefining the origin of time to correspond to 
the time at which the wave swings from leading to trailing. 
Thus without loss of generality we may take k x to be zero. 

When we go on to study the nonlinear development of 
the waves, it is natural to consider disturbances with the 
same coordinate dependence as described above. It is ac- 
cordingly convenient to transform to shearing coordinates 
(x',y'), which are related to (x,y) through 



x = x, y = y + qQ,xt. 



(4) 



The shearing coordinate y specifies the coordinate of an un- 
perturbed fluid element in the direction of the background 
flow, being y at the reference time t = 0. It is also conve- 
nient to work in terms of the velocity perturbation to the 
background shear u — v + qQ,xe y . Then, in shearing coor- 
dinates, the coordinate dependence is on y' alone, so that 
equations Q and Q become 

dp 



dt 

and 
c)u 

Ik 



dy> 



p(u y + qQ,tu x 



du 



= 



+ (u y + qQ,tu x )-^ + 2fl X u — qQ.u x e y 



e y + qQte x \ dP 



dy 



7 + /. 



(6) 



where we have indicated that time derivatives are to be 
taken at constant y' . 



2.2 Lagrangian description 

Equations ([5| and Q describe a form of gas dynamics in 
one spatial dimension. For such problems a simplification 
can often be made by adopting a Lagrangian description 
that uses a spatial coordinate that remains fixed on a fluid 
element. Such a coordinate, yo{y' , t), may be defined through 
the equation 



+ qQtu x 



(7) 



Following Lynden-Bell & Ostriker (19671, we introduce an 



undisturbed 'ghost' flow, which in the present case is simply 
the background shear flow. We consider the mapping of fluid 



elements from the disturbed flow to the ghost flow and take 
yo to be the initial value of y' for the corresponding element 
of the ghost flow. We remark that with this specification we 
do not require that y' for the disturbed flow and yo coincide 
at time t = which in turn allows for a non zero density 
perturbation at t = 0. Note that |7| states that y' moves 
with the component of fluid velocity normal to the phase 
surfaces of constant y' . 

It is convenient to work in terms of the specific volume 
V — 1 1 p. Conversion of the spatial variable from y' to yo is 
effected by using the relation 



dyo 



pa 
P 



V_ 
V ' 



(8) 



which follows from Q and |7f . Here, po is the uniform back- 
ground density, which is fixed for a given fluid element and 
Vo = 1/po- In the Lagrangian description equations Q and 
(|6b thus transform to 



dV 
dt 

du x 

~dt 

8Uy 



d(u y + qQ.tu x ) 
dyo 



dP 



— 2Q,u y = — qCltVa- h f x , 

dyo 

dP 

+ (2 - q)Q.u x = ~Vo^ + fy 
oyo 



(9) 
(10) 

(11) 



2.3 Vortensity conservation 

An equation for the evolution of the ratio of the vertical 
component of vorticity to the density, which we refer to as 
the 'vortensity', is obtained from equations ( |10[ ) and |TTJ by 
multiplying the latter by qQt and subtracting it from the 
former. With the help of equation (|9| we then obtain 



at dyo 
where the vortensity is given by 



Q = (2 - q)QV + Vo 



d(qQ,tu y 



dyo 



(12) 



(13) 



and from now on we omit to indicate that time deriva- 
tives are taken at constant yo- We also find it convenient 
to separate specific volume and the vortensity into back- 
ground and perturbed parts according to V = Vo + SV and 
Q = Qo + SQ, respectively, where Qo — (2 — q)Q/po is the 
background vortensity and 



5Q = (2 - q)Q5V + Vo 



d{q0.tu y — u x ) 
dyo 



(14) 



is the vortensity perturbation. We remark that in the in- 



viscid case, equation (121 shows that vortensity is strictly 
conserved on fluid elements. This remains true when vis- 
cous forces arise from a bulk viscosity, as can be easily seen 
by noting that in this case the viscous forces can be taken 
into account as an addition to the pressure. Furthermore, 



equation (121 is in conservation law form. In Eulerian form 
this reads 



8{pQ) 



dt 



+ 



dy 



]^pQ(u y + qCltu x ) + (f x — qfltfy) 1 ^ = 0. 
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It implies that vortensity is conserved across shocks for any 
kind of small viscosity (even though it may not be con- 
served in passing trough the shock width). Thus for infinites- 
imal viscosity we may adopt SQ — everywhere and apply 
the standard Rankine-Hugoniot conditions to determine the 
changes to other quantities on passing through shocks. Note 
that this is only valid for planar shocks, which we are con- 
sidering here. Curved shocks are in general associated with 
vortensity production (see e.g. Hayesfl957 \ . 

In this context we note that in the linear theory of wave 



excitation described in Heinemann & Papaloizou (2009a 



vortensity perturbations in the form of stationary waves play 
a vital role in that excitation of traveling spiral density waves 
only occurs if such perturbations are present. Wave excita- 
tion happens during a narrow time interval around t — 
(the time of the swing), at which point a pair of counter- 
propagating spiral density waves linearly couples to a sta- 
tionary vortical wave. However, this coupling is not effective 
at later times due to an increasing frequency mismatch be- 
tween the vortical wave and the spiral density waves, so that 
we may take SQ — when describing the dynamics of the 
latter at late times. 



2.4 Conservation of wave momentum in the 
nonlinear regime 

It is also possible to formulate the conservation of wave mo- 
mentum for motions governed by equations |9| to To 
do this we introduce the Lagrangian displacement through 

d$ 



dt 



(15) 



The Lagrangian displacement £ = (£x,£y) is measured rel- 
ative to the background or 'ghost' fluid introduced in Sec- 
tion|2.2|and is not necessarily small. In terms of it, equations 



( 10 \ and (111 are given by 



dt 2 



at 



BP 

2 q n 2 £, x = - q ntv — + f x 



% _ V dP , f 
-dt-- V °Wo +fv ' 



and a time integration of equation |9| yields 
dyo 



sv = Vo- 



(17) 
(18) 



Having expressed the equations of motion in terms of 
Lagrangian displacement, we can now derive a conservation 
law for wave momentum in the form of 



dt dyo dyo ' 

where the wave momentum density is 

dt +QX V dy 
and the associated flux is 
T = -PSV + L. 



U 



(19) 



(20) 



(21) 



Here, C is the Lagrangian density for the system in the ab- 
sence of dissipative forces. This is given by 



E. 



(22) 



where E(V,S) is the internal energy, to which the pressure 
is related by 



P = 



dE 
dV 



(23) 



When there are no dissipative forces, equation ( 19 \ is a strict 



conservation law that implies that the wave momentum flux 
T is constant in a steady state. We remark that this discus- 
sion does not necessarily require a uniform background. It 
applies when Vo is not constant and to the adiabatic case 
with So not constant. By considering the effect of external 
forces the flux T is found to be a momentum flux that con- 
verts to an angular momentum flux on multiplication by an 
assumed radius of the center of the box. 

When dissipation is included wave momentum is no 
longer conserved. From the discussion in Section [23] it is ap- 
parent that we can investigate the role of shocks by including 
a bulk viscosity. This can be done by adding a component 
n to the pressure given by 

n--^ 9 ^ (24) 
v 2 dt ' ( ' 

where £ is the coefficient of bulk viscosity. In this form it is 
clear that this adds a positive contribution to the pressure 
when the gas is being compressed, as is required to support 



shocks. The viscous force corresponding to (24 1 is given by 



/ = -{q£lte x +e y )Vo 



m 

dyo' 



(25) 



When a bulk viscosity is included in this way, equation ( 19 1 
is modified to read 



m_ d{T-nsv) _ n dv 

dt dyo dyo ' 

Thus, wave momentum is no longer conserved. 



(26) 



The right hand side of ( 26 1 is in fact related to the local 



rate of dissipation of energy per unit mass, given by 
,dV 



-n 



dt 



(27) 



Accordingly, the right hand side of equation ( |26[ ) may be re- 
expressed as e(dV/dt)~ 1 dV/dyo- For a non-dispersive travel- 
ing wave, the quantity multiplying e is simply —1/vo, where 
wo is the wave velocity in the shearing coordinate yo- Let 
us consider an outward (inward) propagating trailing wave, 
traveling in the direction of increasing (decreasing) yo, for 
which vo > (vo < 0). As is well known, such a wave carries 
a positive (negative) amount of momentum. Equation ( 26 1 
thus implies that the momentum content of each wave is lost 
to the background. This is consistent with both waves being 
associated with an outward momentum flux. 



3 WEAKLY NONLINEAR THEORY 

Although we are considering sheared disturbances with a 

we have 



2.1l 



specific coordinate dependence (see Section 
so far not made any further approximations. In this section 
we now restrict our consideration to waves with finite but 
sufficiently small amplitude so that all terms of higher than 
quadratic order in the wave amplitude can be neglected. We 
will also generally consider the waves to be in the far trailing 
regime. 
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We start by noting that the number of dynamical equa- 
tions that describe nonlinear spiral density waves can be 
reduced to just a pair of equations by invoking vortensity 
conservation. Indeed, we can dispose of the dynamical equa- 



tion for u y by using ( 141 to express this quantity in terms of 



u x , 5V, and SQ. As we discussed in the end of Section [2. 3[ 
we may set the vortensity perturbation SQ = 0, in which 
case 



1 

qQt 



u x + 



(q-2)Q8- 1 SV 



V 



(28) 



Here, d~o denotes the inverse of the partial derivative with 
respect to yo- An arbitrary function of time arising from 
the integration with respect to yo is left unspecified at this 
point. We will later determine this by enforcing momentum 



conservation (see Section 3.3 below). A consequence of (281 
is that in the limit of late times, the characteristic magnitude 
of u y is smaller than that of u x by a factor of qQt. In this 
limit, it is convenient to adopt the quadratic time variable 



qQt 2 



(29) 



in terms of which the equations of motion governing the 
nonlinear waves are given by 



dV 
dr 

du x 
~dr~ 



= V 



+ 



1 8ut, 



dy qQt dy 



V c 2 dV 2Qu 



+ 



V 2 dyo qQt 



V 



m 

dyo' 



(30) 



(31) 



Here, we recall that recall that c 2 = — V 2 dP/dV evaluated 
at constant entropy is the square of the sound speed and it 
is understood that t = ^/2r/qQ. 

In the limit of large r and vanishing viscosity, equations 



(301 and (311 become equivalent to the standard equations 



for gas dynamics. The terms involving u y lead to small cor- 
rections that arise from the non-uniformity of the shearing 
medium through which the waves propagate. In the linear 
theory given in Hcincma nn fc Papaloizou] ( |2009a| , they are 
seen to induce slow amplitude and phase changes occurring 
on a time scale long compared to the inverse wave frequency. 
These effects appear at first order in a WKBJ treatment and 
we shall treat additional changes arising from nonlinearity 
and viscosity on the same footing below. 

The asymptotic similarity to one-dimensional gas dy- 
namics described in the previous paragraph suggests that 
we describe the nonlinear waves using the Riemann invari- 
ants 



R± = «i± 



V'o 



cdV 



(32) 



in terms of which the equations of motion are given by 



3R± _ ± cVp d 
dr V dyo 



Uy \ 2QUy 



qQt 



Vo 



an 

dyo' 



(33) 



We note that this pair of equations is still exact. Inde- 
pendently of wave amplitude, in the limit of both large r 
and small f, they have solutions consisting of simple waves. 
These are such that the forward propagating wave has 
R+ — and the backward propagating wave has R- = 0. 
The propagation speed cVo/V corresponds to a propagation 
speed c when measured with respect to the Eulerian coordi- 



nate y . In the linear approximation the general solution is 
a superposition of simple waves. 

3.1 An approximate wave equation governing 
weakly nonlinear waves 

At leading order in WKBJ theory, the Riemann invariants 
satisfy first order linear wave equations, according to which 
R± travel in opposite directions at constant amplitude with 
speed Co, being the sound speed c evaluated for V = Vo. 
We consider modifications to the leading order dynamics 
resulting from weak nonlinearity in the wave amplitudes R± , 
corrections resulting from the presence of u y at large r, as 
well as corrections resulting from a small amount of viscosity. 
As the last two corrections are intrinsically small, we only 
consider them to linear order. 

To express quantities of interest in terms of the Rie- 
mann invariants, we note that u x = (R+ + ii_)/2 and 
fy(c/V)dV — (R+ —RJ)/2. From these two relations it 
foil ows that if we set V — Vo + SV, then to first order in 
the wave amplitude 5V we have coSV/Vo ~ (R+ — R-)/2. 
Recalling that Co is the sound speed evaluated for V = Vo, 
the propagation speed is, to the same level of accuracy, given 
by 



cVo 
V 



' Co 



5V\ 
Vo) 



Cl « Co 



R+-R- \ ci 
~2 J co' 



where 



ci 



-V6 



2d(c/V) 



dV 



(34) 



(35) 



v=v 



Note that ci = Co for an isothermal equation of state. 

We focus on the forward propagating wave, which to 
leading order in WKBJ theory is described by i?„ (a cor- 
responding parallel analysis can be done for the backward 
propagating wave). In the dynamical equation for R-, the 
terms containing u y and n (which, as remarked above, we 
consider to be linear) as well as the weakly nonlinear contri- 
bution due to the variable wave speed, given by ( 34 1, involve 
both R- and R + . However, it is possible to argue that R + 
may be neglected in the equation for R-. We can see this 
from the following considerations. 

Assuming that the forward propagating wave is domi- 
nated by R- , we can estimate R+ from its dynamical equa- 
tion, which has small terms proportional to R- as sources. 
In the rest frame of R-, these source terms appear to be 
of high frequency because R+ travels in the opposite direc- 
tions. This leads to the conclusion that the induced R+ will 
be small compared to the original J?_ and so its contribution 
to the small terms in the equation for i?_ may be neglected. 
We note that such a situation is common in weakly nonlin- 
ear hyperbolic systems, see |Hunter fc Keller] ([1983 ) for a 
general treatment. 

With R+ dropped and terms involving u y and n treated 
as linear, the equation for _R_ becomes 

OR- ( ci „ \ dR- 
+ [ C ° + 2T R ~ 



i:)r 



dyo 

1 f 8U V 



qQt V C ° dy 



2Q Ul 



Vo 



m 

dyo'' 



(36) 



where u y and n, respectively defined in ( |28| l and ( |24[ ), arc 
to first order in small quantities given by 
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2qQ.t 



R- + (2 - g)n(c a 2/0 ) _1 R- 



and 

2 Voco or 



(37) 



(38) 



In the last expression for the viscous pressure, we consider 
the bulk viscosity coefficient ( to be a small quantity, in 
which case we may replace dR-/dr by —codR-/dyo- 

The two terms on the right hand side of ( 37 1 give rise to 



various, qualitatively different effects when this expression 
is substituted into Q36n . Both terms lead to slow changes in 
the wave amplitude. The first term also leads to a modifi- 
cation of the propagation speed. The second term also in- 
troduces dispersion. All of these effects are consistent with 
linear WKBJ theory. 

The change in the propagation speed is taken care of 
by transforming to a frame moving with speed 



1 + 



4qQr 



Ct) 



in the yo direction. The new spatial coordinate is thus 



V = Vo 



Co 



1 + 



dr. 



Finally, for notational convenience, we introduce 

ci 



U 



2r, 



R- 



o 



(39) 



(40) 



The governing equation for waves propagating in the yo di- 
rection is then 



dU 
dr 



u 



dU 
dr] 



U_ 
At 



(41) 



where it is understood that d/dr is now to be evaluated 
at constant r\. At this point we remark that as indicated 
above, a parallel analysis can be performed to obtain the 
corresponding equation for the backward propagating wave 
described by U = ciR+/(2co). This is the same as (41 1, 
but with the sign of the dispersive term (first term on the 
right hand side) reversed. Its solutions may be obtained 
from those of (41 1 by use of the transformation r\ — > —rj, 
yo -yo, U -> -U. 

Note that the viscous term increases in magnitude as 
r increases as a consequence of the shear causing the radial 
wavelength of the disturbance to shorten with time. The 
bulk viscosity coefficient £ thus needs to decrease with time if 
the relative importance of the viscous term is not to change. 
The particular form of £ that we use in Section[4]to integrate 
the equation of motion numerically has in fact this property. 



Equation (41 1 may be viewed as a modified form of 



Burgers equation. The first two terms on the right hand 
side correctly account for first order changes in amplitude 
and phase arising in WKBJ theory. The only nonlinearity 
comes from the advection term on the left hand side, which 
causes wave steepening and the formation of shocks. 



discussed the case where nonlinear, dispersive and dissipa- 
tive effects are all small. Physically this corresponds to a 
short wavelength wave in the linear regime with a conserved 
wave momentum which has not had time enough to steepen 
and form shocks. We now go on to consider the case when 
the nonlinear term dominates the dispersive and dissipative 
terms corresponding to a situation where the wave is in the 
process of forming a shock. Setting U = (7/r 1 ^ 4 , adopting 
a new time coordinate f = (4/5) r 5//4 and neglecting disper- 
sive and dissipative terms on the right hand side of (41 1, this 



equation is transformed to the inviscid Burgers equation 



au ^au 

df drj 



0. 



(42) 



It is a well known result in nonlinear acoustics that an 
initially monochromatic wave that satisfies Burgers equation 



steepens to form a sawtooth wave (see Parker 1992 and 
references therein) . Although we do not begin with a strictly 
monochromatic wave we are close to that situation and the 
numerical simulations of Section [4] show that the solutions 
to the full nonlinear equations of motion do in fact approach 
a sawtooth wave, so we shall consider these solutions here. 

A sawtooth wave consists of a periodic array of shocks. 
Because we are considering the inviscid form of Burgers 
equation, we are implicitly taking the limit of vanishing vis- 
cosity. In this limit, the shocks in the sawtooth wave are 
infinitely thin and may be treated as discontinuities. The 
change in U is then determined from the conservation law 



associated with (421. For a stationary (i.e. non-traveling) 
wave, this requires that U changes sign across the shock. 



The solution of ( 42 1 corresponding to a stationary saw- 



tooth wave with wavelength A can be expressed in terms of 
U as 



U = 



5r] 1 

47 1 - (t /t) 5 / 4 



(43) 



where the wave is periodically extended for \rj\ > A/2 and 
ro is an arbitrary integration constant. The wavelength A 
(which we remark is time-independent when measured with 
respect to the shearing coordinate yo) may be specified ar- 
bitrarily (but see below). The solution is linear in rj between 
a periodic array of shocks with amplitude XU/rj. 

Noting that u x = cqU/ci and using the definitions of r 
to express the result in terms of the time t, we see that for 
the fundamental period of the sawtooth wave, when t and r 
are large enough that to may be neglected, we have 



cp 5r/ 
2ci qQ.t 2 ' 



(44) 



The maximum velocity amplitude is obtained by setting 
rj — A/2. The wavelength A thus determines the shock am- 
plitude at a given late time. The fact that this amplitude de- 
cays implies dissipation of energy, which may be interpreted 
as being due to the action of a vanishingly small bulk viscos- 
ity within an infinitely thin shock layer. As we discussed in 
Section |2.4[ dissipation of energy is always associated with 
angular momentum transfer to the background. 



3.2 Nonlinear sawtooth waves 

It is possible to consider solutions of equation ( |41[ ) with dif- 
ferent ordering regimes for the different terms. We have just 



3.3 Correction for dispersion 

We now estimate corrections to the sawtooth solution result- 
ing from dispersion but still take viscosity to be vanishingly 
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small. We work in the limit of large r so that again To in 



(431 is neglected. The governing equation is equation (411 
with C = which reads 



dU 
dr 



dr] At AqftT ^ 

We begin by allowing for a small change to the time 
dependent propagation speed of the wave, which we are re- 
garding to be of the same order as the correction we are esti- 
mating. This gives us the freedom to ensure that the correct 
jump conditions across the shock are maintained. We thus 



add a term C2(r)dU /drj to both sides of (45 1, where C2(r 



is the as yet undefined small increase in the propagation 
speed. In order to deal with the additional term on the left 
hand side, we redefine the co-moving coordinate r) so that it 
is boosted by the speed cz(t) and thus given by 



V = Vo 



Co 



1 + 



AqQr 



dr- 



C2(T 



dr. 



(46) 



In terms of the new co-moving coordinate, equation ( 45 I 
becomes 

^(cod^U , , .8U 
^ " C2(r) a,' 



(W +u dV 
dr dr] 



U_ 
4t~ 



(47) 



where we emphasize again that the last term on the right 
hand side is supposed to lead to a small correction of the 
same order as that produced by the first term on the right 
hand side. 

To find an approximate solution we set U = Uo + U\ 



where Uo is the sawtooth solution given by (431 and Ui is 



a small correction. Dropping terms that are quadratic in Ui 



on the left hand side of (471, and neglecting Ui on the right 



hand side, we find that Ui satisfies the equation 



dUi 
dr 



+ 



57? dUi Ui 
4r drj t 



W + g{r) 



8c r 



(48) 



Here, g(r) is an arbitrary function arising partly from the 
operator d^ 1 and partly from the incorporation of the last 
term on the right hand side of (471. It can be chosen to 



ensure that at any time the wave has no mean momentum, 
i.e. 



A/2 



A/2 



(LT + Ui)d»7 = 0. 



(49) 



To solve (48 1 we set Ui(ri, r) = a(r)r] + /3(t). The functions 



q(t) and f3(r) can be found by substituting this ansatz into 



(481 and integrating with respect to r. In fact, one has only 
to find a(r), since /3(r) can then be determined using the 



constraint (49 1. The dispersive correction is thereby found 
to be 

- 2 (r, 2 - A 2 /12) 



(50) 



16qQrco 

We now turn to the evaluation of C2(r), which we intro- 
duced so as to ensure that the correct shock jump conditions 



are satisfied. From (471 it is apparent that 
U 2 /2 - c 2 U « E#/2 + U {Vi - ca) 

should be conserved when passing through the shock. Since 
Uo changes sign across the shock and U\ does not, this leads 
to 

k 2 X 2 



c 2 (r) = Ui (A/2, r) 



96gf2rco 



(51) 



The corresponding correction to the co-moving coordinate 



defined in ( 46 1 is thus proportional to In r and hence small. 



In order to measure the modification of the sawtooth 
form due to the dispersive correction, we evaluate 

1/2 

r I ' I =- = - (52) 




20 



k 2 A 
15 qQco 



Interestingly, this expression does not depend on time but 
only on the wavelength A. The dependence on A is such that 
the dispersive correction to the sawtooth is small for suffi- 
ciently short wavelengths. However, even for A = 2irco/^l, 
which is the optimal wavelength for excitation by vortensity 
fluctuations in a Keplerian disk ( |Heinemann fc P apaloizou 
|2009a| , we find TZ = tt/15 3/2 w 0.054, suggesting that dis- 
persive corrections to the sawtooth profile may be neglected 
in practice. 



4 NUMERICAL SIMULATIONS 

We have obtained solutions of the set of nonlinear equa- 
tions (JoJ) to (111 by means of numerical simulations. In the 
limit t — > oo this set of equations leads to a pair of equations 
that resembles the standard equations for one dimensional 
inviscid gas dynamics with r = qQ.t 2 /2 as the effective time 
variable. As with those, we expect nonlinear wave steepen- 



ing leading to the formation shock waves (e.g. Landau & 
Lifshitz 1987). In order to represent these, viscous dissipa- 
tion must be included. Although this has to be significant 
only in the thin transition region between the pre-shock and 
post-shock fluid. Here we deal with this by adopting the 



procedure of von Neumann & Richtmyer ( 1950 \ in which an 
artificial viscous pressure, n, is added to the gas pressure, P. 
For the simulations presented here we adopt an isothermal 
equation of state, such that P = c 2 /V, with the isothermal 
sound speed c = Co being constant. Then equations (JoJ) to 
(111 together with (25 become 



po 



dV d(u y + qQ.tu x ) 



at 



dyo 



Po 



— - 2Qu y 



= -qM 



du 



y - + {2- q)Qu a 



d(c 2 /v + n) 

dy 

d(c 2 /v + u) 

dyo 



The form of the viscous pressure n is given by ( 24 1 as 
C dV 



n = 



V 2 dt 



(53) 
(54) 
(55) 

(56) 



the bulk viscosity coefficient to be 
dV 



where, following von Neumann & Richtmyer ( 1950 1, we take 

(57) 



at 



Here, Ayo is the computational grid spacing in yo and ft is 
a constant of order unity which can be adjusted to smear 
shocks over a fixed number of numerical grid points (typi- 
cally between 3 and 5). We adopted /3 = 3 in our simulations. 

We remark that because yo is a shearing coordinate, 
shock widths measured with respect to yo increase linearly 
with time at late times. This can be seen as follows. Inserting 
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(571 into (56 1, and using the fact that for an outgoing wave 



we have to lowest order 

dv n+ dv 

— = -qntc — 
at oy 



(see Section 3.1 1, we find 



n 



(pocPgtttAyof dV 
V dy 



dV 
dyo 



(58) 



(59) 



This expression has an additional factor (qQt) 2 as compared 
to a standard case in non shearing coordinates. Thus the 
shock width scales as PqQ.tIS.yo, which indeed increases lin- 
early with time. Note, however, that yo ~ qiltx at late times, 
which means that the width of a shock remains fixed when 
measured with respect to x. Thus, an observer in the "un- 
sheared" coordinate frame (x, y) will, at late times, see a 
shock of constant width propagating in the x-direction. 

In practice, artificial viscosity is only significant in re- 
gions of strong compression, where dV/dt < 0. Thus, the 
viscous pressure II is often set to zero in regions of rarefac- 
tion, where dV/dt > 0. In our simulations, steep gradients 
never occur in these latter regions, so that this makes lit- 
tle difference, and we have found empirically that slightly 
better results are obtained without this modification. 

Simulations were carried out on an equally spaced com- 
putational grid in yo with typically 16384 grid points. We 
note that because of this rather high resolution, the increase 
of the shock widths as discussed above is not actually visible 
in the plots shown below. The computational domain was 
taken to be periodic with period 2-kc/Q,, which equals the 
optimal wave length for wave excitation for a Keplerian ro- 
tation profile (jHcincmann & Papaloizou 2009a]). In order to 
discretize the equations of motion, we have adopted the same 
staggered leapfrog scheme that was used in | von Neumann~fc| 
|Richtmyer| ( |1950[ ) . This scheme is second order accurate in 
space and time away from shocks, but effectively becomes 
only first order accurate (in space) in their vicinity. 

4.1 Initial conditions 

We start the numerical integration of ( |53[ ) to ( |55[ ) at a time 
corresponding to a couple of orbits before the swing. The ini- 
tial condition consists of a stationary vortical wave of small 
amplitude. To derive this initial condition, we note that in 
the linear inviscid regime, this system of equations is - ac- 
cording to Heinemann & Papaloizou ( 2009a I - equivalent to 



the inhomogeneous second-order wave equations 



dt 2 
and 

&_ 
dt 2 



+ (q 2 Q, 2 t 2 + \)k 22 + k ± 2qQ.ikc 



SV 



T T v 



(ike =F 2fi) p SQ (60) 



+ (q 2 fl 2 t 2 + l)k 2 c + K 



-qQtikyC pot 



(61) 



where SV, 



and the vortensity perturbation SQ are 



assumed to vary harmonically in space as exp(ifcyo)- We note 
that vortensity conservation implies that SQ does not vary 
in time at linear order. 

At early times before the swing, we may drop the double 
time derivative in ( |60[ ) and ( |61[ ) to obtain the slowly varying 
vortical wave solutions 



SV 
V 



(ike =F 2ST2) p SQ 



(q 2 Q. 2 t 2 + l)k 2 c 2 +n 2 ± 2qQikc 



and 

Uy _ 

c 



qQ.tikyC pqSQ 
' (q 2 Q. 2 t 2 + l)k 2 c 2 +k 2 



(62) 



(63) 



where it is understood that the real part of the right hand 
sides is to be taken. 

We initialize the integration with the vortical wave solu- 
tion ( 62 1 and ( 63 1 at t = — 47r/fi, corresponding to two orbits 



before the swing. In this way we generate initial conditions 
that lead to a time dependent evolution that firstly repro- 
duces the linear wave excitation phase described in |Heine-| 
mann & Papaloizou (2009af, in which the vortical wave gives 
rise to the excitation of a pair of counter-propagating spiral 
density waves at t = 0. Subsequent evolution causes the spi- 
ral density waves to enter the nonlinear regime considered 
in this paper and described analytically in Section [3] above. 

We assume a Keplerian rotation profile by setting 
q — 3/2. For the azimuthal wavenumber we choose k = 1/H, 
where H = c/fi is the putative density scale height. This cor- 
responds to an azimuthal wavelength of 2nH, which is the 
optimal wavelength for spiral density wave excitation (see 
Heinemann fe Papaloizou|2009a for details). The amplitude 
of the vortensity perturbation SQ we take to be O.lfi/po- 

4.2 Numerical results 

In Figure[l]we compare the time evolution of solutions to the 
nonlinear equations of motion ( 53 1 to ( 55 1 obtained using 



the artificial viscosity method with solutions of the linear 



wave equations (601 and (61 1. To make this comparison, we 
computed the k — 1 / H Fourier components of the nonlinear 
solutions for u x , u y , and SV (denoted respectively by u x , u y , 
and SV), and plotted the result against the linear solutions. 

The linear and nonlinear solutions are seen to agree 
closely during the leading phase t < and also during the 
early trailing phase t > after the wave excitation has taken 
place. This is consistent with the fact that the amplitude of 
the vortensity perturbation SQ was chosen small enough to 
make the excitation process linear. Nonlinearity occurs later 
as the waves propagate. This is on account of the increas- 
ing amplitude implied as a consequence of wave momentum 
conservation (see Section 2.4 above). A consequence of non- 



linearity and shock development is that the Fourier ampli- 
tudes u x and SV do not continue to increase with time as 
predicted by the linear theory but instead attain maximum 
values and then decay. The maximum relative density per- 
turbation is roughly 20%, which is attained at t ~ 3n/4.Q, 
is similar to that seen in nonlinear simulations of the MRI 
( |Heinemann fc Papa loizou 2009b). It is interesting to note 
that although the nonlinear waves decay while the linear 
waves increase in amplitude they maintain the same phases 
at a given time. We will comment on this further below in 
Section l4~4l 

The decay of the nonlinear wave amplitudes is due to 
a transfer of wave energy to smaller scales due to nonlinear 
mode coupling. That such energy transfer indeed occurs is 
illustrated in Figure|2]where we plot the specific volume V as 
a function of yo, over one wavelength, at three different times 
in the trailing phase. Already half an orbit after excitation, 
nonlinear steepening has lead to a notable distortion of the 
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Figure 1. Evolution of the k = l/H mode amplitudes according 
to the linear (red, attaining the largest values at late times) and 
nonlinear (black) equations of motion. The blue curves represent 
the corresponding vortical wave solutions ( |62[ ) and | |63[ l. The lin- 
ear solution was obtained by integrating the linearized equations 
of motion using a high order Runge-Kutta method. 



initial harmonic profile. After one orbit, further steepening 
has resulted in the formation of two shock fronts correspond- 
ing to the forward and backward propagating waves. These 
separate regions of almost constant specific volume. At this 
time the two shocks travel in opposite directions away from 
the central region. Yet another half of an orbit later, the 
spatial profile remains essentially unchanged, but the shock 
amplitudes have decreased by approximately a factor of two 
because of the energy dissipation associated with them. 

The spatial profiles of the velocity components and the 
specific volume perturbation three orbits after excitation are 
plotted in Figure [3] The shocks in the radial velocity u x are 
seen to be as pronounced as those in the relative specific 
volume perturbation SV/Vo. In this case they separate re- 
gions where the velocity varies almost linearly with yo, cor- 
responding to a sawtooth form. In comparison to SV and 
u x , the spatial variation in the azimuthal velocity field u y 
remains predominantly harmonic. As indicated in Figure [3j 
this sinusoidal variation corresponds to the non-oscillatory 
vortical wave solution that follows from (631. In contrast to 
this, the corresponding amplitudes obtained from (62 1 lead 



to negligibly small values of SV/Vo and u x /c. This suggests 
that in the case of u y , the shock amplitude decays faster than 
the vortical wave amplitude, whereas the opposite seems to 
be the case for SV and u x . This is discussed further below. 




Figure 2. Spatial profile of the specific volume at different times 
in the trailing phase. 
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Figure 3. Spatial profiles of all fluid variables roughly three or- 
bits after excitation (f!t 67r). The black curve corresponds to 
the numerical solution. Also shown in blue are the corresponding 
vortical wave solutions K2[ and 163). 




4.3 Nonlinear damping 

We note that when a forward propagating wave has attained 



a sawtooth form, equation (44 1, with cq — ci — c which ap- 
plies to the numerical simulations, gives u x as a function 
of the co-moving coordinate r\ with a periodic extension for 



r;| > A/2, with A being the azimuthal wave length of the 
original linear wave. The shock discontinuities are accord- 
ingly located at r\ — (n + 1/2) A, where n g Z and the jump 
across the shock is given by 



A(Ua) 



5A 
2qQ.t 2 ' 



(64) 
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Figure 4. Time evolution of the fluctuation energy per unit mass 
£, as defined by equation | |65| l. Black curves not uppermost: En- 
ergy obtained numerically using the artificial viscosity method. 
Red and uppermost curve: Energy obtained from equation l |66| ) 
obtained from weakly nonlinear theory, see Section |4.3| 



This applies to waves traveling in either direction as nei- 
ther is preferred. Thus weakly nonlinear shock jumps decay 
quadratically in time. The decay of the shock amplitudes 
necessarily leads to a decay of the fluctuation energy per 
unit mass, defined as 




_ u 2 x + ul 

2 + VS 



5V 2 
2 ' 



(65) 



Because for weakly nonlinear waves propagating in one di- 
rection, one of the Riemann invariants propagates while the 



other is constant, it follows that (64 1 gives the magnitude 



of the shock jump for both u x and cSV/Vo. From vortensity 
conservation, it follows that u y decays more rapidly than 
SV and u x by one power of t (see discussion in Section [3] 
above) and so may be dropped from (651 as t — > oo. In this 
limit the fluctuation energy averaged over a wavelength as- 
sociated with a pair of shock waves is given by 



{£) = 



1 / 5A 
6 \ 2qQt 2 



(66) 



Thus the fluctuation energy at late time decays as t~ 4 . 

We note that in the limit t — > oo, the above description 
is ultimately independent of the initial wave amplitude. We 
have tested this inherently nonlinear expectation against the 
behavior of our numerical simulations. In Figure [4] we plot 
the fluctuation energy ( 65 I for waves initiated with four dif- 



ferent amplitudes of the vortensity perturbation 5Q. Accord- 
ing to the linear excitation mechanism described by |HeintT| 
mann & Papaloizou (2009a I, this amplitude, determines the 



amplitude of the excited linear wave. For the examples illus- 
trated the fluctuation energies obtained shortly after excita- 
tion differ by almost two orders of magnitude. However, as t 
approaches large values, the fluctuation energies associated 
with the numerical solutions all approach the predicted late- 



time fluctuation energy (66 1. This is because the more the 



initially linear wave is amplified, the sooner it becomes non- 
linear, and larger amplitude waves decay faster than smaller 
amplitude ones. 



4.4 Shock propagation 

In this section we will investigate the space-time dependence 
of the nonlinear waves. We remind the reader that for the 
forward propagating weakly nonlinear solution described in 



Section 3.1 the co-moving coordinate is given by (46 1 as 
V = yo-cr- — JJ 1 + , (67) 



AqQ. 



where we have now reverted back to quadratic time vari- 
able r used in Section |3.1| For the backward propagating 
wave c is replaced by — c in (67 1. Apart from asymptotically 



insignificant amplitude (through A) and logarithmic correc- 
tions, the functional form of the shock wave solution is the 
product of a time-dependent amplitude and a function of 
the phase j/o — cr as is also the case for linear waves. Thus 
the speed at which the shock fronts propagate is approxi- 
mately equal to the sound speed, being the phase speed of 
linear waves in the far trailing regime. 

We are now in a position to understand the evolution of 
individual Fourier harmonics as obtained from the numeri- 
cal solution of the nonlinear equations, as shown in Figure]!] 
There we see that while nonlinear effects do lead to a de- 
crease in the oscillation amplitude, they do not seem to alter 
the oscillation period. This follows immediately from the fact 
that the linear and nonlinear waves maintain approximately 
the same phase as a function of time. 



4.5 Shock wave profile 



In Section |4.2| we illustrated the spatial profile of a pair of 
shock waves observed in a numerical simulation. This is char- 
acterized by a top-hat profile in the specific volume pertur- 
bation and a double sawtooth profile in the radial velocity. A 
similar double sawtooth profile is also seen in the azimuthal 
velocity, but there it is disguised by the spatially sinusoidal 
vorticity perturbation that we used as initial condition. As 
we will now demonstrate, all of these features can be under- 
stood from weakly nonlinear theory. 

In Figure|5]we illustrate schematically the superposition 
of two counter-propagating sawtooth waves, corresponding 
to the weakly nonlinear asymptotic solutions derived in Sec- 
tion [3] The arrows indicated the direction of propagation. 
Note that the forward and backward propagating wave have 
the same slope. As we noted earlier, this is because the latter 
is obtained from the former by simultaneously reversing the 
sign of the wave amplitude and the spatial coordinate. Since 
the two sawtooth waves have the same slope, subtracting one 
from the other yields a top hat profile for the specific volume 
perturbation. Conversely, adding the two waves results in a 
double sawtooth profile for the radial velocity perturbation. 

We can gain additional confidence in the weakly nonlin- 
ear theory by testing it directly against the numerical data. 
In Figure [6] we plot analytic solutions on top of numerical 
solutions at a time corresponding to roughly three orbits af- 
ter excitation. To achieve the best possible agreement, we 
have incorporated in the analytic solution both the disper- 



sive correction to the sawtooth wave given by ( 50 \ as well 



as the stationary vortical wave that we used as initial con- 
dition. 

As Figure [6] shows, the analytic solutions are essen- 
tially indistinguishable from the numerical ones. We note 
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Figure 5. Sketch of the superposition of two counter-propagating 
sawtooth waves. The arrows indicated the direction of propaga- 
tion. Note that the waves travel towards a region of low pressure 
(since 8p ~ 1/8V). 



V/V 





i 


i i 


1.02 








1.01 








1.00 








0.99 


i 


i i 






u x /c 




0.02 


i . 


i i 


0.01 






0.00 








0.01 






0.02 


i i 


" 1 






Uy/c 





0.004 - 
0.002 - 
0.000 ' 
-0.002 - 
-0.004 - 



-A/2 A/2 

yo 

Figure 6. Comparison of wave forms roughly three orbits af- 
ter excitation (fit s« 67r). Black: Numerical solution. Dotted red: 
Weakly nonlinear solution. 

that the only free parameter that enters the analytic solu- 
tions is the location of the shock front. The magnitude of the 
jump across the shock is in fact correctly predicted by ( |64[ ). 
We also note that a superposition principle evidently holds. 
This confirms what we have argued for in Section |3.1| i.e. 
that weakly nonlinear, counter-propagating sawtooth waves 
do not interact to a good approximation. 



Close inspection of the wave profile of the specific vol- 
ume perturbation between shock fronts reveals that it is 
slightly curved. This curvature is a dispersive effect that is 
correctly captured by the dispersive correction to the leading 
order sawtooth solution. Unlike for SV, there is hardly any 
curvature in the wave profile of the radial velocity u x . This 
is expected, as it turns out that the dispersive correction, 
which is quadratic in the spatial coordinate, enters R+ and 
R- with opposite sign, so that it cancels out when their sum 
is taken to obtain u x . We note that these effects cannot be 
attributed to the vortical wave as its amplitude is less than 
1% of the shock amplitude. In addition, the curvature in the 
specific volume perturbation profile is everywhere concave. 
The vortical wave on the other hand would give rise to to a 
sinusoidal variation that is concave in one half of the domain 
and convex in the other. 

Also illustrated in Figure [6] is the wave profile of the 
azimuthal velocity u y . The analytic solution shown consists 
of a superposition of the vortical wave and the weakly non- 
linear shock waves. The wave form is mostly sinusoidal, in- 
dicating that the vortical wave dominates. This is because 
in the case of u y , the vortical wave decays as t~ , whereas 
the shock waves decay as t -3 , so that the former will always 
dominate at late times. 

For SV and u x , the vortical wave and the shock waves 
decay equally fast (as t~ 2 ). In general, one might therefore 
expect to see the vortical wave at late times. Since this is 
not the case, we conclude that shortly after the initial wave 
excitation, the shock wave should be much larger in ampli- 
tude than the vortical wave. Inspection of Figure [I] shows 
that this is indeed the case. 



5 DISCUSSION AND APPLICATION TO MRI 
SIMULATIONS 

We have considered the propagation of density waves in a 
shearing box. These were assumed to be functions only of 
the shearing coordinate y + qQxt and time t. The waves 
were presumed to have been excited by vortensity fluctua- 
tions produced by MRI turbulence as described by |Heine-| 
|mann fc Papaloizou| ( 2009a b). In this process it is only the 
form of the vortensity fluctuation at the time when waves 
swing from being leading to trailing that is significant. Just 
after the swing inward and outward propagating waves are 
assumed, as found in simulations, to be in the linear regime. 
Subsequently, as the waves propagate, their wavelengths in 
the x direction decreases, while the conservation of wave 
momentum results in increased amplitude, causing them to 
enter the nonlinear regime and the formation of shocks. 

In Section[3]of this paper we developed an analytic the- 
ory that is applicable to the weakly nonlinear regime. We 
derived a modified Burgers equation governing the dynam- 
ics of weakly nonlinear waves, which contains terms describ- 
ing nonlinear steepening, dispersion, and viscous diffusion. 
We obtained nonlinear sawtooth solutions with weak shocks 
valid for late times and estimated corrections resulting from 
dispersion. 

In Section [4] we presented numerical solutions of the 
nonlinear equations without a small amplitude approxima- 
tion. These led to sawtooth waves with a profile that was in 
excellent agreement with that obtained from the weakly non- 
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Figure 7. Pseudo-color images of the vertically integrated mass 
and momentum densities taken from the simulation described in 
|Heinemann fc Papaloizou| l |2009b| |. 



linear theory at late times, in particular with regard to rate 
of decay of the shock velocity jump being ultimately pro- 
portional to t~ 2 . The solutions constructed from the theory 
were composed of a combination of forward and backward 
propagating waves, resulting in a double sawtooth profile for 
the radial velocity u x and a top hat profile for the specific 
volume perturbation SV. 



5.1 MRI simulations 

An important issue is the extent to which the description of 
weakly nonlinear waves we have provided above applies to 
the density waves excited in MRI simulations. These waves 
are found to be ubiquitous in such simulations. As a typi- 
cal example of what is found in many realizations we show 
in Figure [jj a snapshot of the (vertically integrated) mass 
and momentum densities found in a simulation described 
in Heinemann & Papaloizou (2009b I. From this figure it is 



apparent that the wave structure is significantly more pro- 
nounced in u x and 5p/po ~ — SV/Vo than in the azimuthal 



velocity 



This is expected from our discussion in Sec- 



tionary vortical wave. In addition, Figure [jj indicates a top 
hat like profile in SV . To examine this more closely, Fig- 
ure[8]shows plots of the fluid variables taken at a fixed value 
of x. These may be compared to the corresponding plots 
in Figure [5] It is seen that although there are superposed 
fluctuations in the MRI simulation, there is good agreement 
between the profiles: double sawtooth for u x and top hat 
for 8V (corresponding to an inverted top hat for the profile 
of Sp). This correspondence is less clear for u y , but this is 
not unreasonable on account of the expected dominance by 
vortical perturbations. 

It should be noted that in an MRI simulation, spiral 
density waves propagate through a field of turbulent fluc- 
tuations, and it might be expected that the waves inter- 
act with these fluctuations. This could in general lead to 
wave front distortion and energy loses resulting from random 
secondary density wave generation through interaction of 
the primary wave with turbulent eddies (eg. |Lighthill|1953 




-0.10 



Figure 8. Cut through Figure [jj along y at a fixed x RJ 0.15 (no 
average over x). 



Howe|| 1971a|b||Ffowcs Williams fc Howe|l973"l ) . We also re 



mark that a spiral density wave is expected to undergo wave 
front distortion as it propagates through a medium with lo- 
cal turbulent velocity and/or local sound speed fluctuations 
(eg. Hesselink fc Sturtevant|1988 1. This may occur in a sta- 
tionary medium with no generation of secondary waves with 
different frequency. 

We comment that a description of losses through wave- 
turbulence interactions by an effective turbulent viscosity is 
rather problematic. Such a description necessitates statisti- 
cal averaging over many realizations of the system, see e.g. 
Balbus & Papaloizou ( 1999k. Thus it most likely describes 



the possible variation that different realizations can display 
rather than the behavior of any one realization, where dis- 
tortion of the wave profile rather than diffusion may occur, 
an issue that leads to interpretational difficulties, see |Ffowcs| 
|Williams & Howe ( 1973 ) for an extensive discussion. These 



authors also point out that an effective turbulent viscos- 
ity description fails to meet the reasonable expectation that 
shock waves should be influenced mostly by small scale tur- 
bulent motions. 

In support of the above view, we comment that ex- 



perimental results presented by Plotkin & George ( 1972 \ 



for wave profiles associated with sonic booms propagating 
through a turbulent atmosphere indicate good agreement 
with non turbulent theory apart from random perturbations 
to the wave profile and a thickening of the shock front (that 



remains narrow on a macroscopic scale) (sec also Giddings 
et al.|[2001 l. Although the context differs, this view would 



appear to be consistent with the results presented in Fig- 
ures [7j and 

Nonetheless, it is likely that interaction with disorga- 
nized turbulence causes the waves to decay more rapidly 
than predicted by our weakly nonlinear wave theory. An in- 
spection of the results in |Heinemann &: Papaloizou| ([2009b ) 
indeed indicates a decay rate faster than t~ . However, this 
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may be affected by the particular choices of Reynolds and 
Prandtl numbers that are made for numerical tractability. 
In this context it is important to note that the decay rate 
produced in the simple weakly nonlinear wave theory enun- 
ciated here gives, what we expect from the discussion above, 
to be a lower limit as far as the simulations are concerned 
and it does not differ greatly from what is actually observed. 
Thus the description of excited waves given by the simple 
theory and the current simulations is unlikely to be greatly 
modified when the transport coefficients are changed. 
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